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Abstract 

In this paper we present a one dimensional second order accurate method to solve Elliptic equations 
with discontinuous coefficients on an arbitrary interface. Second order accuracy for the first derivative is 
obtained as well. The method is based on the Ghost Fluid Method, making use of ghost points on which 
the value is defined by suitable interface conditions. The multi-domain formulation is adopted, where 
the problem is split in two sub-problems and interface conditions will be enforced to close the problem. 
Interface conditions are relaxed together with the internal equations (following the approach proposed 
in [TU] in the case of smooth coefficients), leading to an iterative method on all the set of grid values 
(inside points and ghost points). A multigrid approach with a suitable definition of the restriction oper- 
ator is provided. The restriction of the defect is performed separately for both sub-problems, providing 
a convergence factor close to the one measured in the case of smooth coefficient and independent on the 
magnitude of the jump in the coefficient. Numerical tests will confirm the second order accuracy. 
Although the method is proposed in one dimension, the extension in higher dimension is currently un- 
derway |12j and it will be carried out by combining the discretization of |10) with the multigrid approach 
of 11 for Elliptic problems with non-eliminated boundary conditions in arbitrary domain. 

Introduction 

Elliptic equations with jumping coefficients across a one-codimensional interface T arise in several applica- 
tions. Let us mention as examples the steady-state diffusion problem in two materials with different diffusion 
coefficient separated by an arbitrary interface, the Poisson equation coming from the projection method in 
incompressible Navier-Stokes equation for fluids with different density, the porous-media equation to model 
the oil reservoir, electrostatic problems, and many others. In order to close the problem, interface conditions 
related to the jump of the solution and of the flux across the interface are included. In all these problems 
the interface may be arbitrary (not aligned with a line grid) and can change in time. 

Numerous techniques have been developed to treat such problem. Interface-fitted grid methods such as 
the ones based on Finite Elements Methods O [5] are not suitable in case of moving interface, because a 
re-meshing grid is needed at each time step and this makes the computation expensive. Then an approach 
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treating the interface embedded in a Cartesian grid and moving according to the velocity field of the fluid is 
preferred. Since the interface may not be aligned with the grid, a special treatment is needed. The simplest 
method makes use of the Shortley-Weller discretization [30J, that discretizes the Laplacian operator with 
usual central difference away from the interface, and makes use of a non symmetric stencil in the points close 
to the interface, adding extra-grid points on T. While jumping condition on the solution is straightforward 
to discretize on interface points, the jump in the flux (involving the normal derivative) cannot be immediatly 
discretized in more than one dimension. In fact, Shortley-Weller discretization requires that the value of the 
normal derivative of the solution on both sides of the interface is suitably reconstructed at the intersection 
between the grid and the interface. This approach is adopted, for example, by Hackbusch in |20j to first 
order accuracy, and by other authors (see j6] and references therein) to second order accuracy. However, 
the method proposed by Bramble in [B] for second order accuracy is quite involved and not recommendable 
for practical purposes. 

Methods based on embedding the domain in a Cartesian grid without adding extra-grid points are 
derived from the pioneering work of Peskin [26], where the Immersed Boundary Methods is introduced to 
model blood flows in the heart. In that paper a source term is localized on the the boundary and the 
method makes use of a discretized delta-function, leading to a first order accuracy. A second order accurate 
extension to jump coefficients is the Immersed Interface Methods, first developed by LeVeque and Li in [23]. 
Such method uses a six-point stencil to discretize the elliptic equation in grid points close to the interface T 
and the coefficients of such stencil are found by Taylor expansion of the solution. Jump conditions on the 
interface are then used to modify the coefficients appearing in the equation corresponding to nodes near T, 
in such a way that the overall discretization is second order accurate. Non-homogeneous jump conditions 
are allowed on the function and on the normal flux. 

Another method which achieves second order accuracy by modifying standard difference formulas was 
proposed by Mayo in [24] for solving Poisson or biharmonic equation on irregular domains. Such method 
embeds the irregular domain in a regular region with a Cartesian grid and discretizes the equation on the 
whole region, by suitable extension of the solution outside. 

In all these methods the only unknowns are the values on the grid points and the stencil may cross the 
interface, leading to a quite involved procedure to reach the desired accuracy, since the derivative of the 
solution may jump crossing the interface and values from the other side are used in the computation. 

A rather simple method to use standard five-point stencil even close to the interface is the Ghost-Fluid 
Method, introduced by Fedkiw et al. in [14J. Here the authors point out that a two-phase problem could 
be reduced in two sub-problems by a multi-domain formulation, and each sub-problem may be discretized 
with the same technique used to solve a single problem with Dirichlet/Neumann boundary conditions. Such 
method makes use of extra grid points (ghost points) outside the domain in order to keep unchanged the 
symmetry of the stencil even for inside points close to the interface. In ghost points, interface conditions 
are enforced in order to close the discrete system. 

Methods based on ghost points are discussed in [16], where Gibou et al. proposed a second-order 
accurate method for Dirichlet conditions on regular Cartesian grid. The value at the ghost nodes is assigned 
by linear extrapolation, and the whole discretization leads to a symmetric linear system, easily solved by a 
preconditioned conjugate gradient method. A fourth order accurate method is also proposed in |17j . Other 
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methods use a non-regular Cartesian grid, such as in [9] , where Gibou et al. present finite difference schemes 
for solving the variable coefficient Poisson equation and heat equation on irregular domains with Dirichlet 
boundary conditions, using adaptive Cartesian grids. One efficient discretization based on cut-cell method 
to solve more general Robin conditions is proposed by Gibou et al. in [25], which provides second order 
accuracy for the Poisson and heat equation and first order accuracy for Stefan type problems. 

Other approaches based on cut-cell methods obtained by a Finite Volume discretization are presented 
in |22j . Cells that are cut by the boundary requires a special treatment, such as cell- merging and rotated-cell, 
in order to avoid a too strict restriction of the time step dictated by the CFL condition. 

Several methods have been also proposed to model the interaction between multiphase flows and solid ob- 
stacles, such as Arbitrary Lagrangian Eulerian (ALE) |15tll3j. Distribute Lagrangian Multiplier (DLM) |18j . 
penalization methods [29 , 2J. In [8j a combination of penalization and level-set methods is presented to solve 
inverse or shape optimization problems on uniform Cartesian meshes. In [32] Zhou et al. proposed a Matched 
Interface and Boundary (MIB) method for elliptic problems with sharp-edged interfaces. 

In time-dependent problems requiring the solution of an elliptic problem at each time step an itera- 
tive solver is preferred with respect to a direct problem, since a good initial guess (the solution at the 
previous time step) is provided. Most iterative method for jumping coefficient are based on Domain Decom- 
position Methods [27], either with or without overlapping. Such methods are based on the multi-domain 
formulation, i.e., the problem is split in two sub-problems and interface conditions are enforced to achieve 
two sub-problems with respectively Dirichlet and Neumann boundary (coupled) conditions on the inter- 
face/boundary. Each sub-problem is solved and the solution at the interface is used to provide an updated 
right-hand side for the other sub-problem, and so on iteratively. A drawback of this method is that associa- 
tion between the Dirichlet/Neumann boundary condition and the sub-domain cannot be arbitrary (see [27\ 
pag. 12]). 

Most applications require second order accuracy in the gradient: for example, in projection method for 
incompressible Navier-Stokes equation, the gradient of the pressure is used to correct the fictitious velocity 
field leading it to satisfy the free-divergence condition. Also high-order accuracy [T7] may be required, 
for instance when turbulence and shock interact, or high frequency wave propagation are presented in 
inhomogeneous media [1]. 

In [10] a second-order accurate discretization for elliptic problems in arbitrary domain and mixed bound- 
ary condition is provided, together with a convergence proof for the iterative solver for first order accuracy. 
The method is based on transforming the stationary problem into a fictitious evolutionary problem, both 
inside the domain and on the boundary. The problem is then discretized on a regular grid using non elim- 
inated boundary conditions to determine the proper relaxation equation for the ghost points. The whole 
procedure is made efficient by a multigrid technique, as illustrated in 

The present paper provides a second order discretization of the problem based on the ghost-point method 
on regular Cartesian grid described in [10] and makes use of an iterative solver whose convergence is speeded 
up by a multigrid approach [11]. Interface conditions are neither eliminated from the discrete system (they 
are strongly coupled and their elimination is too hard to perform in more than one dimension) nor directly 
enforced (which leads to a non-convergent iterative method): they are relaxed together with the interior 
equations. This leads us to an iterative scheme for the set of all unknowns (internal points and ghost points). 
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The method works also for non- homogeneous interface conditions. Although this paper provides a ID 
description of the method, the generalization of the approach in higher dimension is currently underway [12] 
and can be obtained in an almost straightforward manner combining results from |XQ^ 111) . 

Several multigrid approaches exist in literature to treat the jumping coefficient problem in 2D when 
the interface is aligned with the Cartesian grid. We mention the method based on operator-dependent 
interpolation [TJ [21] , where the interpolation is carried out by exploiting the continuity of the flux instead 
of the gradient of the solution, and the method based on Galerkin Coarse Grid Operator [28J, which makes 
the algebraic problem more expensive from a computational point of view and does not take advantage from 
the fact that the discrete problem comes from a continuous problem. 

In our approach we use the standard interpolation operator and discretize the operator in the coarser 
grid in the same way as in the fine grid, without making use of Galerkin conditions. But, since the defect 
may jump crossing the interface, a separated restriction for both sub-problems is needed, as performed 
in for arbitrary domain with mixed boundary condition (without jumping coefficient). This approach 
provides a good convergence factor, comparable with ones measured for no-jumping case. We also show 
that the convergence factor does not depend on the magnitude of the jump in the coefficient. Interface 
conditions are relaxed, then have to be transferred to the coarse grid as well. In one-dimensional case this 
task is trivial, since such conditions are just two real values that can be copied to the coarse grid. In 
higher dimension interface conditions are stored in ghost points, which can show a complex structure for 
arbitrary interface. The restriction of interface condition defect can be carried out in the same manner of 
the restriction of boundary condition defect described in [11] for problems with non-eliminated boundary 
conditions: the defect is first extrapolated outside the domain and then transferred to the coarse grid in the 
same manner as the restriction of the defect of inside equations, i.e., without using values from the other 
side of the boundary. This work is currently underway. 

The rest of the paper is divided in 3 sections. In the first section we describe the second order accurate 
discretization of the model problem and the iterative scheme obtained by the relaxation of the interface 
conditions. The second section is devoted to the multigrid approach, with a care description of the transfer 
operators. In section 3 some numerical test is performed, to show the second order accuracy in the solution 
and in its first derivative as well. We measure also the convergence factor and compare it with the convergence 
factor obtained by other methods. 

1 Second order accurate discretization 

In this section we obtain a second order accurate numerical method to solve an elliptic equation with 
discontinuous coefficients. After introducing the model problem, we provide a discretization and an iterative 
solver of the linear system. In some applications one may be interested in second order accuracy also for the 
derivative of the solution. In numerical tests of Sec. [3] we show that the method is second order accurate in 
the solution and in its first derivative. 
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1.1 Model problem 

Let us consider the model problem 



dx 
u(0) 



[0,1], 



(1) 



U(l) = fifl. 



where the diffusion coefficient 7: [0, 1] — >• M jumps on an interface a €]0, 1[, i.e., is a smooth function in 
[0, a[ and in ]a, 1], but may be discontinuous across a. We assume 7 > e > in all the domain. If we solve 
this problem by standard central differences on a uniform grid, the accuracy of the method degrades to first 
order. 
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Fig. 1: Computational domain Q with an arbitrary interface a. 



Let 
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U =U\[ 0>a [, U" = u\] a> i], 7 =7l[Q,a[ ; 7"=7l]a,l] 

be the restriction functions of the solution and of the coefficient on the two subdomains. We split the 



problem into the following subproblems: 
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u R (l) 


= 9i- 





(2) 



(3) 



In order to close the problem, we must provide an additional boundary condition for each of u and u R on 
the interface a. This additional conditions are inferred to the requirement that the solution u and the flux 
ju' are continuous across a. Introducing the jumping operator on a 



w 



lim w 

x— >a+ 



lim w, 



the additional boundary conditions may be resumed as 

[u] = 0, [yu'] 
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and are called transmission conditions [27]. They can be inferred by a physical requirement: for instance, in 
steady-state diffusion problems in two materials, the temperature and its flux are required to be continuous 
across a. Non-homogeneous interface conditions may appear, for example, in presence of a delta-function on 
the right hand side / = f\ + S a , with f± E C°([0, 1]). Precisely, the two following problems are equivalent: 

u(0) = g , u(l)=gi, 

~{y^)=fi'm [0,a[U]a,l] 
ax ax 

u(0)=g , u {l)=gi 
[u} = 0, [ 1 u']=-C. 

In the following we suppose the right-hand side is a regular function in the two sub-regions, and non- 
homogeneous interface conditions are allowed: 



u 



9D, [iu'] = g N . (4) 



Such general case is relevant for some applications, for example pressure equation for incompressible flow in 
presence of surface tension at the interface. 

The two subproblems ^ and ^ are then coupled on a and cannot be solved separately. The whole problem 
becomes 

'' 'l Rd} t) = /in I'M] W 



dx \ dx 

u L (0) = g , u R (l) = 91 (7) 
9D, [ju') = 9N- (8) 



u 



1.2 Discretization 

Let ./V be an integer, h = 1/{N + 1) be the spatial step and xo, x%, . . . , xn, x^r+i be the equally spaced grid 
points, with Xj = j h. Let J be such that xj < a < xj+i (see Fig. [I]). We write J = [a[, where [ - L 
denotes the integer part. We will denote by Lj[w] the quadratic interpolant of w in nodes {xj-\,Xj,Xj+x}. 
By Uj [u R ] we denote the component of the numerical solution which approximates u L (xj) [u R (xj)], while 
we intend fj = f(xj), = -f L (xj), jf = j R (xj). 

Let us discretize the system ([7]). Discretizing Eq. ^ on nodes xi, x%, . . . , xj using central differences for 
the solution u L and linear interpolation for the coefficient function 7 , we obtain: 

h bh W - u ^ + 7 ivi ( u f - u J+i)) = u> j = 1, ■ ■ • j, (9) 
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where 7^ +1 / 2 = (7^ + 7^+1 )/2- l n Eq. (^) for j = 1 the value Uq is given by the Dirichlet condition (|7j): 
Uq = go- It can be easily eliminated from ([9|, but we will leave it in the system just for simplicity. The 
same applies for discretizing Eq. Q in node xjy. 

Eq. ([9]) for j = J needs to know the values of u L and 7 L in node xj + \. Since v! and 7 are discontinuous, we 
cannot use respectively u R +l and 7j +1 , because this may result in a loss of accuracy, since it smears out the 
coefficient 7 and the numerical solution itself, while both jump on the interface. Then we need to add an 
additional grid point value for the numerical solution u L (x j+i), called ghost point value, and to extrapolate 
j L up to the first ghost point Xj+x- The same argument holds for u R and 7^ in their ghost point xj, when 
discretizing Eq. ^ in node xj + \. 

The unknowns of the numerical method are therefore the N + 4 quantities 

u%,...,u5 +1 ,u$,...,u% +1 . (10) 

This approach has been called Ghost Fluid Method and used in the context of multi-fluid flows [H] . The two 
additional unknowns Uj +1 and u R require two additional boundary conditions to close the system, which are 
given by the transmission conditions Q, resulting in a 2 x 2 sub-system. We will not solve this sub-system 
for Uj +l and u R , but we instead leave it in the whole linear system, which will be solved iteratively. The 
extrapolation for the coefficient functions j L and 7^ is simple linear extrapolation: 

75+x = 2 7^ - 7^1 , 7 j = 2 7j+i - 7j +2 • 

Using then central differences to discretize ([5]) and (pj), linear and quadratic interpolation to discretize 
respectively the two conditions ([8]), we obtain the following second order (N + 4) x (N + 4) linear system: 

4 = 90 (11) 
1 (jf_ h (uf - uU) + («f - uf +l )) = Si J = 1, ..- J (12) 

((1 _ # )U R + $ U R +l ) _ ((1 _ )U L + = g D (13) 

7* Lfj[u R ](a) - 7 £ Lj-i [«*](«) = g N (14) 

p («f " uf-i) + 7f + | («f - uf +1 )) = j = J + 1, . . . N (15) 

t$ + l = <7i, (16) 

with 7^ and 7^ obtained by linear interpolation: 



-y£ = (1 - 0)7* + 7* = (1 " #hj + #7* 



J+l 



and?? = {a-xj)/he [0,1]. 

If we apply a simple iterative method such as Gauss-Seidel or Jacobi to this linear system, in general it 
will not converge, unless we solve the 2x2 sub-system of transmission conditions, eliminating them from 
the whole system. This elimination is easy to perform in one dimension, but becomes quite involved in 
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higher dimension. Therefore, we prefer to work with the whole linear system without eliminate transmission 
conditions from it, in order to extend the method to higher dimension in a forthcoming paper |12j . Then 
we have to find a different approach to solve iteratively the previous linear system. This can be done by 
relaxing the transmission conditions. 



1.3 Iterative method 



In order to find a convergent iterative method to solve the linear system ( 11 )-{ 16 ) , following the approach 
introduced in [10] we solve the associate time- dependent problem in the unknowns u L (x,t) and u R (x,t) for 
(x,t) G [0,1] x (0,+oo): 



u L {0,t) 
du L 



ot 



9o 



du L 



Ot 

du R 



dx \ dx 

du 

1- 



x G [0, a[ 



dt 



dx 

t*D {9d - [«]) 



9n 



du R 
~dT 



d 



u R (l,t) 



^ dx \ dx 
91- 



(17) 
(18) 

(19) 

(20) 

(21) 
(22) 



where fi is a positive function, and [id and are two positive constants, that will be set in Sec. |1.4| to 
satisfy some stability condition. 

The choice of the sign of the two constants \xr> and /^tv is crucial and requires some explanation. Roughly 
speaking, when replacing a vector equation F(w) = for F : IR m — > W 1 by u)= F(u) u, we have to be sure 
that the solution is asymptotically stable, i.e. that A(V^F) < 0. Eq. (20) will be used to compute u 



H 
J ' 

therefore the derivative of the right hand side of Eq. (20) with respect to Uj has to be negative, to ensure 



convergence to equilibrium. Eq. (19) is used to determine Uj +1 by a transport equation on u (x,t). Since 



xj + \ > a the propagation speed hn1 L associated to u L (x,t), has to be positive. 
We are obviously interested in the steady-state solution and the time t represents an iterative parameter. 



We observe that transmission conditions (19) and (20) can be replaced by 



du R 



dt 
du L 



Mat 9n 



7 



du 
dx 



dt 



A»r>([«] ~9d) 



because both choices lead to the same steady state conditions. 

To obtain a second order accurate solution in space we are allowed to discretize first order accurate the 
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time derivative. Using forward Euler in time and central differences in space for (18) and (21), we obtain 
(superscripts L and R are omitted): 



(m+l) 



(m) 



7i_i «■ 



(m) (m) 



+ 7 



i+i 



11 



(m) (m) 
U j+l, 



h 2 



(23) 



where j = 1, . . . , J for u L and j = J + 1, . . . , N for u R . Choosing the maximum time step allowed by the 
CFL condition for diffusion equation, i.e., fij At = /i 2 /(7j+i/2 + 7j-i/2)> Eq. (23) becomes: 

' ( t l2 i ( m ) i ( m )\ 



It 



(m+l) 



(24) 



'3-2 ^ 'j+2 

where j = 1, . . . , J for -u 1 ' and j = J + 1, . . . , N for u R . Observe that such equation is the one obtained by 
applying Jacobi iteration to Eqs. (12) and (15). 

Let us discretize Eq. (19). The time derivative is discretized by forward Euler at the ghost point xj+i, 
which is the quantity we want to compute. The jump is discretized as in (14), so it is second order accurate. 
We obtain the iteration: 



u 



L,(m+1) 

J+i 



^ + ^Ai ( 7 * h'j[u R ^}(a) - ^ L'j^lu^ia) - g N ) 



Likewise, in Eq. (20) we discretize the time derivative in xj, obtaining: 



u 



R,(m+1) 
J 



R,(m) 



+ (j,d At ((1 - 0) u y m) + *uj+ x ) L >™ - (1 - $)uf (m) + 0u?ff> + !1D 



(25) 



(26) 



Iterations (24), (25) and (26) constitute the iterative scheme to solve problem ([!]) to second order accuracy. 
1.4 Choosing constants \xb and fi^ for transmission conditions 



In (25) and (26) two arbitrary constants \id and hn appear. Following the same argument as in [TU], such 
constants will be chosen in order to satisfy some stability condition for the equation where they appear. 
This procedure is not rigorous because it does not take into account the coupling between the equations, 
and does not consist in a convergence proof. However, in all numerical tests we performed, the conditions 
we find seem to guarantee convergence. 

Constant pLjj is introduced in Eq. (20), which is just a relaxation of the jump condition. Then we require: 

Hd At < 1. (27) 

This condition will ensure positivity, and is a factor 2 more stringent than just stability restriction. For 
practical purpose, we set \xd At = 0.9. In order to obtain a condition on fijy, we rewrite Eq. (19) as follows 
(we have supposed for simplicity homogeneous jump = 0): 



~dt 



+ Mat 7 



du L 
dx 



mi 



dx 



t G (0,oo). 
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This is a simple convection equation with speed /j.nJ L '■ Then a simple CFL condition for convection equation 
might be 

h 

UN At < —jr. 

Numerical experiments show that this condition is not enough, especially in the case ^ R /^ L ;§> 1. An 



explanation of this behavior may be that the right-hand side of (28) is not stationary when the convection 



evolves in time, but it depends on time itself by u R . An acceptable condition is 

» N At< / (29) 

max {7,7} 

For practical purpose we choose fi^ At = 0.9 hj max {7^,7^}. Numerical tests show that conditions (27) 



and (29) are sufficient for guarantee convergence, but not necessary. A more detailed analysis is in progress. 
Notice that fiAt = 0(h 2 ), fi^ At = 0(h), fj,£> At = 0(1). Furthermore, only the product of the constants 
times At enters into the conditions, therefore we may imagine that At = 1. 



2 Multigrid approach 



The convergence of the iterative method proposed in Sec. 1.3 is usually very slow. To accelerate the 
convergence we use a multigrid strategy. To make the iteration scheme (24)-(26) a building block for 



an efficient multi-grid solver, we must be sure that such iteration (relaxation scheme) has the smoothing 
property, i.e. that after few steps, the error becomes smooth (not necessarily small). Roughly speaking, 
the high-frequency components of the error reduce quickly. We do not explain all multigrid features, but 
just what is different from classical multigrid approach, remanding to the literature for more details (e.g., 



see [3T1 [T9l [T] ) . The iteration scheme (24)-(26) is a Jacobi-like scheme, as mentioned in Sec. 1.3 Jacobi 
scheme is not a good smoother, since high-frequency components of the error reduce slowly. A good smoother 



is instead the Gauss-Seidel scheme. Then, we use a Gauss-Seidel version of (24)-(26) as relaxation scheme, 
i.e. 

1 

L,(m) 



L,(m+1) 



ft* 



L,(m+T) 



L,(m) 



j = l,...,J 
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u 



L,(m+1) 
J+l 

fl,(m+l) 
J 



uy^ + m At[^h'j[u R ^)( 



< > 



i[u L ](a) 



u 



R,(m) 



R,(m+1) 



updated in t 



+ ^ D At (^(1 - 
1 



$) Uj 



L,(m+1) 



+ du 



L,(m+1) 
J+l 



(1 



\ R,(m) 
V)Uj 



where in (|31|) we intend u L such that 



re same order reported in (10) 



2 . R,(m+1) . R,(m) 



u L,(m+l) f or j < J _)_ 1 an( 4 u L 



9N 



a RAm) 

+ 9D 



j = J+l,.. 



(30) 
(31) 

(32) 
(33) 



J+i 



L,(m) 
L J+1 ■ 



The unknowns are 



10 



In order to explain the multigrid approach, we just describe the two-grid correction scheme (TGCS), because 
all the other schemes, such as V-cycle, VK-cycle, F-cycle or Full Multigrid cycle, can be easily derived from 
it (see [3H Sections 2.4 and 2.6] for more details). Let us introduce some notation. For a grid of spatial step 
h, we denote: 



J 



a 
h 



J 



S(Qh) = { w h = (w L ,w R ) such that w L : {xo, . . . , xj+i} — > M, w R : {xj, . . . , xn+i} — > ^} 
S (Qh) = { w /i = (w L , w R ) such that w L : {x\, . . . , xj} ->■ R, w R : {xj+i, . . . , xn} M} 

u h = ((Uj)j=0,...,J+l, (u R )j=J,...,N+l) G S(Qh) 

lh = {{ij)i=Q,...,j+i,(rif)j=j,...,N+i) e S(Q h ) 

o 

fft G5 such that fft(zj) = /j 
L ft : 5(f2 h ) x 5(O ft ) — >S (^h) such that 



(L fc (7 fc ,Ufc)) j 
(L fc (7 ft ,Ufc)),- 



l)+Tf + jW-«i+i)) J 



^ 2 V 7f_ i («f " uf_ x ) + 7 f+i («? - uf +1 ) ) if J > J + 1 



\u ■■ s(n h ) 



such that 

l-D /7i .QN..H , „q„„R \ [(■> „q\„.L , „□„,£ 



[n h ]h = ((1 " W + faf+i) ~ ((1 " W + M+i) 
[ • , ■ ]^ : S(fi ft ) x 5(n h ) — »• R such that 
hh, Ufc]? = 7 R L'j[u R }(a) - ^ L'j^[u L ](a) 



The linear system (!!)-( 16) can be resumed as follows: 



L h (j h ,u h ) 
[ u h\ h 



f/, 

Si- 



(34) 
(35) 
(36) 
(37) 
(38) 



For simplicity we assume that N + l = 1/h is a power of 2. The TGCS consists into the following algorithm: 
1. Set initial guess u/j = 0. 
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2. Relax v\ times on the finest grid: for k from 1 to v\ do (30), (31), (32) 

o 

3. Compute the defects €S {^h), 9d,9n £ ^ : 



*h = fft + L h (j h ,u h ) 
9D 



3d - [u h ] h 



lA' 



4. Transfer the defect to a coarser grid with spatial step 2/i by a suitable restriction operator 

*2h = 4* ( r />) • 

5. Solve exactly the residual problem on the coarser grid in the unknow e 2 h £ S(£l2h) 

Lh(j2h,e2h) = *2h 

[e2h\ h = 9d 

n2h,U2h\h = 9N 



4 








c (7V+l)/2 

6. Transfer the error to the finest grid by a suitable interpolation operator 

e h = ll h {e 2h ). 



7. Correct the fine-grid approximation 



u h = u h + e h . 



8. Relax u 2 times on the finest grid: for k from 1 to v 2 do (30), (31), (32) 



To complete the description of TGCS, we have just to explain the steps concerning grid migration (steps 4 
and 6). 

2.1 Transfer grid operators 

In this section, we describe the transfer grid operators for vertex-centered grid. Observe that coefficients j L 
and j R can be transferred in an exact manner by a simple injection operator. 
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2.1.1 Restriction operator 



T R 

Since such operator will act on the defect = (r^,r^) ES {^h) (step 4), we perform the restriction from 
a fine grid to a coarser grid separately for r^ and r^. This is justified by the fact that the defect r^ of the 
left domain may be very different (after few relaxations) from the defect r^ of the right domain, especially 
in the case of high jumping coefficient, i.e., max {7^/7^, 7^/7^} >> 1. In addition, these defects are very 
different also from the defects of jumping conditions go and gw, because the operators scale with different 
power of h. 

Let us describe the restriction of r^ by the operator (-^h)^ ( see Fig. [2). Let xj be the closest grid point to 
a from the left in the fine grid (see Fig. [TJ. Let x be a grid point of the coarse grid. If x < xj we will use 
the standard full- weighting restriction operator (FW): 



l2k) L r L h (x) 



(r h (x -h) L + 2 r h (x) L + r h (x + h) L ) 



while if x = xj we reduce to an upwind linear convex combination from the left direction: 

L 



= "I r h( x ) + (1 " wi)r£(* - h), 



(39) 



(40) 



since in x + h only r^ is defined and not rf. In our tests we found that uj% = 1/2 gives better results than 
vi = 3/4. 

The operator (i^J^ works in a similar manner: let xj+i the closest grid point to a from the right in the 
fine grid. If x > xj+i we will use the standard full- weighting restriction operator (FW): 



1 



r h (x -h) R + 2 r h (x) H + r h (x + h) 



R 



while if x = xj we reduce to an Upwind mean value from the left direction: 

1 



The whole restriction reads 



4^ 



(41) 



(42) 



jh\ L r L (jh\ R T R 



In the upper part of Fig. [2] is represented the case in which we have to use (40) and (41). The only other 
possible case is that we have to use (39) and (42). 



2.1.2 Interpolation operator 

Since such operator will act on the correction &2h = ( e 2h> e 2h) e S(&2h) (step 4), we perform the interpolation 
from a coarse grid to a finer grid separately for e^ and e^ (see middle and lower part of Fig. [2]), but always 
using the standard linear interpolation: 
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Restriction operator 



Q 



h 



t — t — * — t — i — t- 



Q 



i — + — k — + — i — + — t 



a 



Interpolation operator 



A 1 h 




Q 



2h 



Interpolation operator 



H f 



A — + — k — + — * — + — *r 



a 



1 



Fig. 2: Fine and coarse grid for transfer operators. The dashed lines represent the action of the restriction (top) and the 
interpolation (middle and bottom) operators. 



i^f^h^j) = e 2h( x j) if j is even 

( J f) L< 4(*i) = H e 2h( x J-i) + e 2h( x 3+i)) if i is odd. 



( 7 f ) Re 2h( x j) = e 2h( x j) if j is even 

{ll h ) R 4h( x j) = |(e&(x i - 1 ) + ^(x i+1 )) ifjisodd. 



2^o„, - ( ( T?h\ L L ( T 2h\ R -R 



The whole interpolation reads 

Remark. 1 (Coarser operator) We observe that the discrete operator L<ih on the coarser grid (step 
5) is just the operator obtained discretizing directly the continuous operator in the grid with spatial step 
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2h, and not the operator obtained by Galerkin condition 

L 2h = I2I1 L h Ih 



(45) 



The last approach, typical of algebraic multigrid, makes the algebraic problem more expensive from a 
computational point of view and does not take advantage of the fact that the discrete problem comes from 
a continuous problem. 

Remark. 2 (y-cycle) The F-cycle algorithm is easily obtained from the TGCS recursively, namely 
applying the same algorithm to solve the residual equation in step 5. To terminate the recursion, an exact 
solver is used to solve the residual problem when the grid achieves a fixed level of coarsening. We denote 
by V(ui, z^-cycle the y-cycle performed with v\ pre-relaxations and V2 post-relaxations. 

Remark. 3 (W-cycle) The VF-cycle is similar to the T^-cycle, with the only difference that the residual 
problem is solved recursively two times instead of one (in general schemes, 5 times, but 5 > 2 is considered 
useless for practical purposes). 



3 Numerical tests 



In this section we confirm numerically the second order accuracy of the discretization of Sec. 1.2 and compute 
the convergence factor p of the multigrid approach for several examples, to confirm the independence of p 
from the spatial step h and the magnitude of the jumping coefficient. 

Second order accuracy is gained also for first derivative of the solution, as it is shown by the comparison 
between exact first derivative and the numerical derivative obtained by central difference of the numerical 
solution. 

In all numerical tests, we choose an arbitrary interface a G]0, 1[ and an analytical expression of the exact 



solution u = (u ,u H ) and of diffusion coefficient 7 



Then we reconstruct the data /, gr> and , 



perform the multigrid technique, and compare the numerical solution with the exact solution to compute 
the order of accuracy by the slope of the best-fit line. In all our tests we use the following stopping criterion 
for the V— cycle 



u 



(m+l) 



U 



(m) 



U 



(m+l) 



< TOL. 



This will ensure that the actual relative error satisfies 



(m+l) 



< 



\ e h\ 



TOL 
'~p 



The tolerance we used is TOL = 10 -6 , which ensures that the error in the solution of the algebraic system 
is always lower than truncation error. For each example we show a table in which we list the errors, and 
the value in the third [fifth] column and z'-th row of the table indicates the accuracy order, computed as 
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log 2 (ei-i/ei), where is the L°°-error of the numerical solution [derivative] indicated in the second [fourth] 
column and i-th row. 

To compute the asymptotic convergence factor, we use the following estimate: 



P = P 



(m) 



„(ro-l) 



which is reliable for m large. In order to avoid difficulties related to numerical instability due to machine 
precision, we will always use the homogeneous model problem as a test when we want to compute the 
asymptotic convergence factor, namely Eq. ([!]) with / = go = g\ = and homogeneous jump conditions, 
and perform the multigrid algorithm starting from an initial guess different from zero. Since in this case 
we are just interested in the convergence factor and not in the numerical solution itself (which approaches 
zero), a reasonable stop criterion will be 

\ n ( m ) — rt ( m - 1)1 

& — A — ' < io- 2 . 

Several tests are performed for each example, based on the different size of the finest and coarsest grids. 
The finest grid is obtained dividing the domain [0, 1] into N + 1 intervals, while the coarsest grid is obtained 
dividing the domain into N c + 1 intervals. 



3.1 Example 1 

We choose (see Fig. [3]) 



a = 0.343, 



= sin(57nr) 
~2 



3 + cos(5"7rx) 
10 9 (10 + sin(57rx)) 



Fig. [4] shows the numerical results and the second order slope of the best-fit line for the L°°-error of the 
numerical solution and its derivative. Table [I] shows the convergence factor for different values of N and N c . 



3.1 



Table 1: Measured V(l, l)-cycle convergence factor for the numerical test of Ex. 
number of grid points in the finest grid; N c + 2 number of grid points in the coarsest grid. 



We use N + 2 



N+l 

Nc+1 


32 


64 


128 


256 


512 


1024 


2048 


4096 


16 


0.15 


0.16 


0.15 


0.17 


0.19 


0.15 


0.15 


0.15 


32 




0.14 


0.12 


0.17 


0.19 


0.15 


0.15 


0.15 


64 






0.07 


0.16 


0.19 


0.15 


0.15 


0.15 


128 








0.11 


0.17 


0.15 


0.15 


0.15 
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O exact solution 
3 r • numerical solution 



2.5 -<% 

© 







© 



1.5 - ® 



1 3 



© 

® 

© 
© 




25 

20^2 
15 



10 
5 

-5 
-10 
-15 
-20 
-25 



O exact first derivative 
. numerical first derivative 



© 



© ® 
© 

© 
©® 
© 




15r 



x 10 



* ydiffusion coefficient 



10 




0.5 



x 10 



10 



-10 



-15 L 



* source term f 




1 



0.5 



Fig. 3: We refer to Ex. 3.1 The data are computed for N = 64. 



Error in the solution 
bestfit slope=-2.00 



0) 

o 




■f^ Error in the first derivative 
bestfit slope=-2.00 

OF 




N + l 


u — 




oo 


order 


||u'- 


<l 


oo 


order 


64 


1.87 


•10- 


-2 




3.33 


•10- 


-1 




128 


4.59 


•10' 


-3 


2.03 


8.38 


•10- 


-2 


1.99 


256 


1.13 


•10" 


-3 


2.02 


2.12 


•10- 


-2 


1.98 


512 


2.77 


•10- 


-4 


2.03 


5.35 


•10- 


-3 


1.99 


1024 


6.95 


•10- 


-5 


2.00 


1.34 


•10- 


-3 


2.00 


2048 


1.73 


•10- 


-5 


2.01 


3.36 


•10- 


-4 


1.99 


4096 


4.48 


•10- 


-6 


1.95 


8.21 


■10" 


-5 


2.03 


8192 


1.11 


•10- 


-6 


2.01 


2.07 


•10- 


-5 


1.99 



Fig. 4: We refer to Ex 



3.1 



Left: Representation of the L x 
derivative. The slope of the best-fit lines is respectively s = 
errors and order of accuracy computed by subsequent errors. 



-error of the numerical solution and its 
-2.00 and s = -2.00. Right: List of 



3.2 Example 2 

We choose (see Fig. [5]) 



a = 0.743, 



= sin(57nr) 
~2 



, R 



3 + cos(5-7rx) 
10 9 (10 + sin(57rx)) 
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The only difference with respect to the previous example is the value of a. 

Fig. [6] shows the numerical results and the second order slope of the best-fit line for the L°°-error of the 
numerical solution and its derivative. Table [2] shows the convergence factor for different values of N and N c . 



O exact solution 
. numerical solution 



5 
1 3 



© 

© 

© 
© 




© © 



© 
© 



© 
© 



O exact first derivative 
. numerical first derivative 




10 



15 r 



# ydiffusion coefficient 



10 



; 



0.5 



10 



10 



* source term f 



-10 



-15 L 



0.5 



1 



Fig. 5: We refer to Ex. 



3.2 



The data are computed for N = 64. 



Error in the solution 
bestfit slope=-2.00 



°-4 

CD 
O 




Error in the first derivative 
■ bestfit slope=-2.00 

OF 




JV + 1 


u — 




oo 


order 


||u'- 


U 'J 


oo 


order 


64 


1.86 


■10- 


-2 




3.38 


•10' 


-1 




128 


4.63 


•10" 


-3 


2.01 


8.38 


■10" 


-2 


2.01 


256 


1.15 


•10" 


-3 


2.01 


2.10 


•10" 


-2 


2.00 


512 


2.86 


•10" 


-4 


2.01 


5.26 


•10- 


-3 


2.00 


1024 


7.24 


•10" 


-5 


1.98 


1.30 


•10- 


-3 


2.01 


2048 


1.80 


•10" 


-5 


2.01 


3.28 


•10- 


-4 


1.99 


4096 


4.48 


•10" 


-6 


2.01 


8.21 


•10- 


-5 


2.00 


8192 


1.12 


•10" 


-6 


2.00 


2.05 


•10" 


-5 


2.00 



3.2 



Left: Representation of the L° 



Fig. 6: We refer to Ex 
derivative. The slope of the best-fit lines is respectively s = 
errors and order of accuracy computed by subsequent errors. 



error of the numerical solution and its 
-2.00 and s = -2.00. Right: List of 
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Table 2: Measured V(l, l)-cycle convergence factor for the numerical test of Ex. 3.2 We use N + 2 
number of grid points in the finest grid; N c + 2 number of grid points in the coarsest grid. 



N+l 

Nc+1 


32 


64 


128 


256 


512 


1024 


2048 


4096 


16 


0.13 


0.13 


0.11 


0.14 


0.15 


0.15 


0.15 


0.15 


32 




0.13 


0.11 


0.15 


0.15 


0.15 


0.15 


0.15 


64 






0.11 


0.13 


0.15 


0.15 


0.15 


0.15 


128 








0.15 


0.15 


0.15 


0.15 


0.15 



3.3 Example 3 

We choose (see Fig. [7]) 



a = 0.283 



u 



,R 



= sin(57rx) 
~2 



10 9 (10 + sin(5vrx)) 
3 + cos(57rx) 



Fig. [8] shows the numerical results and the second order slope of the best-fit line for the L°°-error of the 
numerical solution and its derivative. Table [3] shows the convergence factor for different values of N and N c . 



O exact solution 
. numerical solution 




25 — 
20 £> 
15 
10 
5 



O exact first derivative 
• numerical first derivative 



© 



-10 - 



© 



© 

© 

© ® 
© 



© 

15 - O© 

® 



X 10 

15, 

* -/diffusion coefficient 



10 



0.5 



10r 



x 10 



J2 



I 

* source term f 



# 
5 _% 



0** 



0.5 



Fig. 7: We refer to Ex. 



3.3 



The data are computed for iV = 64. 



3.4 Example 4 

We choose (see Fig. [9]) 



, u L = e x2 / 7 L = 10 9 (10 + sin(5vrx)) 

a - U.6L6, <j uR = eSin{5nx) , I 7 R = 3 + cos ( 57rx ) 
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y^r Error in the solution 
bestfit slope=-2.01 




^ Error in the first derivative 
bestfit slope=-2.00 

OF 




N + l 


u — 




oo 


order 


||u'- 


<\ 


oo 


order 


64 


2.07 


■10" 


-2 




3.15 


•10" 


-1 




128 


5.15 


•10" 


-3 


2.01 


7.76 


•10" 


-2 


2.02 


256 


1.20 


•10" 


-3 


2.10 


1.90 


•10- 


-2 


2.03 


512 


2.19 


•10" 


-4 


2.46 


5.57 


•10" 


-3 


1.77 


1024 


6.10 


•10" 


-5 


1.84 


1.33 


•10" 


-3 


2.07 


2048 


1.76 


•10" 


-5 


1.79 


3.09 


•10" 


-4 


2.11 


4096 


5.05 


•10" 


-6 


1.80 


7.60 


•10" 


-5 


2.02 


8192 


1.22 


•10" 


-6 


2.05 


1.86 


•10" 


-5 


2.03 



3.3 



Left: Representation of the L° 



Fig. 8: We refer to Ex 
derivative. The slope of the best-fit lines is respectively s = 
errors and order of accuracy computed by subsequent errors. 



-error of the numerical solution and its 
= -2.01 and s = -2.00. Right: List of 



Table 3: Measured V(l, l)-cycle convergence factor for the numerical test of Ex. 



3.3 



We use N + 2 



number of grid points in the finest grid; N c + 2 number of grid points in the coarsest grid. 



N+l 

Nc+1 


32 


64 


128 


256 


512 


1024 


2048 


4096 


16 


0.09 


0.10 


0.12 


0.15 


0.15 


0.15 


0.15 


0.15 


32 




0.09 


0.10 


0.15 


0.15 


0.15 


0.15 


0.15 


64 






0.12 


0.15 


0.15 


0.15 


0.15 


0.15 


128 








0.13 


0.15 


0.15 


0.15 


0.15 



Fig. 10 shows the numerical results and the second order slope of the best-fit line for the L°°-error of the 
numerical solution and its derivative. Table [4] shows the convergence factor for different values of N and N c . 



3.5 Independence of convergence factor from the jump in the coefficient 

In this section we show that the convergence factor does not depend on the jump in the coefficient. We 
choose 



a 



0.543, 



u L = f 7 L = 10 p 



ii 



R 



' 



7 



R 



1 



and start the multigrid process with an initial guess different from zero, in order to compute the asymptotic 
convergence factor. We list the results in Table [5} 
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O exact solution 
• numerical solution 




^ Error in the solution 
— bestfit slope=-1.99 




O exact first derivative 
■ numerical first derivative 



x 10 

15| 

* ydiffusion coefficient 
10 



x 10 



10 




© 



0.5 




Fig. 9: We refer to Ex. 3.4 The data are computed for N = 64. 



^ Error in the first derivative 
— bestfit slope=-2. 00 




JV + 1 


u — 




oo 


order 


||u'- 


<l 


oo 


order 


64 


1.59 


•10- 


-2 




3.46 


•10" 


-1 




128 


3.99 


•10" 


-3 


2.00 


8.74 


•10" 


-2 


1.98 


256 


9.66 


■10" 


-4 


2.05 


2.22 


•10- 


-2 


1.98 


512 


2.23 


•10" 


-4 


2.12 


5.74 


•10- 


-3 


1.95 


1024 


5.25 


•10" 


-5 


2.08 


1.47 


•10- 


-3 


1.97 


2048 


1.68 


•10" 


-5 


1.64 


3.31 


•10- 


-4 


2.15 


4096 


4.12 


•10" 


-6 


2.03 


8.36 


•10- 


-5 


1.98 


8192 


9.87 


•10- 


-7 


2.06 


2.13 


•10" 


-5 


1.97 



Fig. 10: We refer to Ex 



3.4 



Left: Representation of the L°°-error of the numerical solution and its 
derivative. The slope of the best-fit lines is respectively s = —1.99 and s = —2.00. Right: List of 
errors and order of accuracy computed by subsequent errors. 



Remark. (Comparison with Domain Decomposition Method) 

Domain Decomposition Method (DDM) is another iterative method to solve elliptic problems with discon- 
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Table 4: Measured V(l, l)-cycle convergence factor for the numerical test of Ex. 3.4 We use N + 2 
number of grid points in the finest grid; N c + 2 number of grid points in the coarsest grid. 



N+l 

Nc+1 


32 


64 


128 


256 


512 


1024 


2048 


4096 


16 


0.17 


0.12 


0.14 


0.18 


0.17 


0.15 


0.16 


0.15 


32 




0.11 


0.14 


0.16 


0.15 


0.15 


0.15 


0.15 


64 






0.06 


0.14 


0.15 


0.15 


0.15 


0.15 


128 








0.12 


0.15 


0.15 


0.15 


0.15 



Table 5: Measured V(l, 1) asymptotic convergence factors for a problem with a jumping coefficient 
of the order 1(F 



p 





1 


2 


3 


4 


5 


p 


0.11 


0.10 


0.11 


0.11 


0.11 


0.10 



tinuous coefficient, based on solving iteratively the two subproblems 



d_ ( L du L '( m+1 ) 
dx \ ' dx 

u L,(m+l)( ) 

u L,(m+l)( a ) 



,d_ ( R du R '( m+1 ) 
dx \ ' dx 



in [0, a[ 



9o 

u R ^ m \a) 



7 



R du R '( m + 1 '> (a) 



f m]a,l] 



7 
91 



dx 



(46) 



(47) 



until convergence. A little drawback of this method is that, in order to guarantee the convergence, it must 
be a > 0.5 (see |27} pag. 12]). Our method may be regarded as a DDM, but in place of solving a subproblem 
to provide the right-hand side for the other subproblem (and so on iteratively), we just perform a relaxation 
on a subproblem, and with the guess obtained we build the right-hand side of the other subproblem, as it 
can be seen in Sec. 



1.3 With this relaxing strategy, the convergence is always guaranteed, as showed in 



numerical tests. 



Conclusion 

A second order discretization for elliptic equation with discontinuous coefficient on an arbitrary interface 
has been provided. Second order accuracy in the derivative is obtained as well. The linear system is solved 
by an iterative method obtained relaxing the interface conditions. The iterative method is then speeded 
up by a proper multigrid approach, which transfers separately the defect for both sub-problems obtained 
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from the multi-domain formulation. The measured convergence factor is close to the one measured in the 
case of smooth coefficients and it does not depend on the magnitude of the jump in the coefficient. The 
method is similar to Domain Decomposition Methods, but a single relaxation sweep is performed in each 
subdomain instead to solve it completely. This makes the method more flexible and there is no restriction 
on the relative size of the two subdomains. 

This paper is the building-block for a future work in higher dimension [12] , which will be carried out by 
combining the second order discretization in arbitrary domain with smooth coefficients [10] and the multi- 
grid treatment of problems with non-eliminated boundary conditions in arbitrary domain . 
Other future works concern the convection-diffusion equation in a moving domain, in order to study appli- 
cations modeled by a Stefan- Type problem. A level-set function will keep track of the moving interface. 
All this extensions will be coupled with the use of Adaptive Mesh Refinement to obtain accurate solution 
in the case of domain with complex boundary. A proper multigrid approach is under investigation for all 
these works. 
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